home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 7 / Apprentice-Release7.iso / Environments / PowerLisp 2.01 / PowerLisp 2.01 ƒ / Library / random.lisp < prev    next >
Encoding:
Text File  |  1996-05-22  |  6.8 KB  |  216 lines  |  [TEXT/ROSA]
  1. ;;; -*- Mode: Lisp; Package: Kernel -*-
  2. ;;;
  3. ;;; **********************************************************************
  4. ;;; This code was written as part of the CMU Common Lisp project at
  5. ;;; Carnegie Mellon University, and has been placed in the public domain.
  6. ;;;
  7. ;;;
  8. ;;; **********************************************************************
  9. ;;;
  10. ;;; Functions to random number functions for Spice Lisp
  11. ;;;
  12. ;;; Originally written by David Adam.  Python tuning, better large integer
  13. ;;; randomness and efficient IEEE float support by Rob MacLachlan.
  14. ;;;
  15. (in-package "COMMON-LISP")
  16. (provide :random)
  17. (export '(random-state random-state-p random *random-state*
  18.       make-random-state))
  19.  
  20. (in-package "POWERLISP")
  21. (export '(%random-single-float %random-double-float random-chunk
  22.                    random-fixnum-max))
  23.  
  24. ;;;; Random state hackery:
  25.  
  26. (defconstant random-const-a 8373)
  27. (defconstant random-const-c 101010101)
  28. (defconstant random-max 54)
  29. (defconstant single-float-digits (float-precision 1.0))
  30. (defconstant double-float-digits (float-precision 1.0))
  31. (defconstant word-bits 32)
  32. (defconstant single-float-significand-byte (byte 52 0))
  33. (defconstant double-float-significand-byte (byte 52 0))
  34.  
  35. ;;; Inclusive upper bound on the size of fixnum kept in the state (and returned
  36. ;;; by random-chunk.)  Must be even.
  37. ;;;
  38. (defconstant random-upper-bound (- most-positive-fixnum 3))
  39. (defconstant random-chunk-length (integer-length random-upper-bound))
  40. (deftype random-chunk () `(integer 0 ,random-upper-bound))
  41.  
  42. (defun fixnump (x) (typep x 'fixnum))    ;; PowerLisp addition
  43.  
  44. (defvar rand-seed 0)
  45. (defstruct (random-state
  46.         (:constructor make-random-object)
  47.        ;;#| (:make-load-form-fun :just-dump-it-normally) |#
  48.         )
  49.   (j 24 :type index)
  50.   (k 0 :type index)
  51.   (seed (make-array (1+ random-max) :initial-contents
  52.             (do ((list-rands () (cons (rand1) list-rands))
  53.              (i 0 (1+ i)))
  54.             ((> i random-max) list-rands)
  55.               (declare (fixnum i))))
  56.     :type simple-vector))
  57.  
  58.  
  59. ;;; Generates a random number from rand-seed.
  60. (defun rand1 ()
  61.   (declare (optimize (inhibit-warnings 3)))
  62.   (setq rand-seed
  63.     (mod (+ (* rand-seed random-const-a) random-const-c)
  64.          (1+ random-upper-bound))))
  65.  
  66.  
  67. (defvar *random-state* (make-random-object))
  68.  
  69.  
  70. (defun copy-state (cur-state)
  71.   (let ((state (make-random-object
  72.         :seed (make-array 55)
  73.         :j (random-state-j cur-state)
  74.         :k (random-state-k cur-state))))
  75.     (do ((i 0 (1+ i)))
  76.     ((= i 55) state)
  77.       (declare (fixnum i))
  78.       (setf (aref (random-state-seed  state) i)
  79.         (aref (random-state-seed cur-state) i)))))
  80.  
  81. (defun make-random-state (&optional state)
  82.   "Make a random state object.  If State is not supplied, return a copy
  83.   of the default random state.  If State is a random state, then return a
  84.   copy of it.  If state is T then return a random state generated from
  85.   the universal time."
  86.   (cond ((not state) (copy-state *random-state*))
  87.     ((random-state-p state) (copy-state state))
  88.     ((eq state t) (setq rand-seed (get-universal-time))
  89.               (make-random-object))
  90.     (t (error "Argument is not a RANDOM-STATE, T or NIL: ~S" state))))
  91.  
  92. ;;;; Random entries:
  93.  
  94. (declaim (special start-block random %random-single-float %random-double-float
  95.               random-chunk))
  96.  
  97. ;;; random-chunk  --  Internal
  98. ;;;
  99. ;;; This function generates fixnums between 0 and random-upper-bound, 
  100. ;;; inclusive.  For the algorithm to work random-upper-bound must be an 
  101. ;;; even positive fixnum.  State is the random state to use.
  102. ;;;
  103. (declaim (ftype (function (random-state) random-chunk) random-chunk))
  104. (defun random-chunk (state)
  105.   (let* ((seed (random-state-seed state))
  106.      (j (random-state-j state))
  107.      (k (random-state-k state))
  108.      (a (- (- random-upper-bound
  109.           (the random-chunk
  110.                (svref seed
  111.                   (setf (random-state-j state)
  112.                     (if (= j 0) random-max (1- j))))))
  113.            (the random-chunk
  114.             (svref seed
  115.                (setf (random-state-k state)
  116.                  (if (= k 0) random-max (1- k))))))))
  117.     (declare (fixnum a))
  118.     (setf (svref seed k)
  119.       (the random-chunk (if (minusp a) (- a) (- random-upper-bound a))))))
  120.  
  121. ;;; %RANDOM-SINGLE-FLOAT, %RANDOM-DOUBLE-FLOAT  --  Interface
  122. ;;;
  123. ;;;    Handle the single or double float case of RANDOM.  We generate a float
  124. ;;; between 0.0 and 1.0 by clobbering the significand of 1.0 with random bits,
  125. ;;; then subtracting 1.0.  This hides the fact that we have a hidden bit.
  126. ;;;
  127. (declaim (inline %random-single-float %random-double-float))
  128. (defun %random-single-float (arg state)
  129.   (declare (type (single-float (0f0)) arg)
  130.        (type (or random-state null) state))
  131.   (* arg
  132.      (- (%make-single-float
  133.      (dpb (ash (random-chunk (or state *random-state*))
  134.            (- single-float-digits random-chunk-length))
  135.           single-float-significand-byte
  136.           (%single-float-bits 1.0)))
  137.     1.0)))
  138. ;;;
  139. (defun %random-double-float (arg state)
  140.   (declare (type (double-float (0d0)) arg)
  141.        (type (or random-state null) state))
  142.   (let ((state (or state *random-state*)))
  143.     (* arg
  144.        (- (%make-double-float
  145.        (dpb (ash (random-chunk state)
  146.              (- double-float-digits random-chunk-length
  147.             word-bits))
  148.         double-float-significand-byte
  149.         (%double-float-bits 1d0))
  150.        (logxor (ash (random-chunk state)
  151.             (- word-bits random-chunk-length))
  152.            (ash (random-chunk state)
  153.             (- random-chunk-length word-bits))))
  154.       1d0))))
  155.  
  156. ;;;; Random integers:
  157.  
  158. ;;; Amount we overlap chunks by when building a large integer to make up for
  159. ;;; the loss of randomness in the low bits.
  160. ;;;
  161. (defconstant random-integer-overlap 3)
  162.  
  163. ;;; Extra bits of randomness that we generate before taking the value MOD the
  164. ;;; limit, to avoid loss of randomness near the limit.
  165. ;;;
  166. (defconstant random-integer-extra-bits 10)
  167.  
  168. ;;; Largest fixnum we can compute from one chunk of bits.
  169. ;;;
  170. (defconstant random-fixnum-max
  171.   (1- (ash 1 (- random-chunk-length random-integer-extra-bits))))
  172.  
  173.  
  174. ;;; %RANDOM-INTEGER  --  Internal
  175. ;;;
  176. (defun %random-integer (arg state)
  177.   (declare (type (integer 1) arg) (type random-state state))
  178.   (let ((shift (- random-chunk-length random-integer-overlap)))
  179.     (do ((bits (random-chunk state)
  180.            (logxor (ash bits shift) (random-chunk state)))
  181.      (count (+ (integer-length arg)
  182.            (- random-integer-extra-bits shift))
  183.         (- count shift)))
  184.     ((minusp count)
  185.      (rem bits arg))
  186.       (declare (fixnum count)))))
  187.  
  188. (defun random (arg &optional (state *random-state*))
  189.   "Generate a uniformly distributed pseudo-random number between zero
  190.   and Arg.  State, if supplied, is the random state to use."
  191.   (declare (inline %random-single-float %random-double-float))
  192.   (cond
  193.    ((and (fixnump arg) (<= arg random-fixnum-max))
  194.     (rem (random-chunk state) arg))
  195.    ((typep arg 'single-float)
  196.     (%random-single-float arg state))
  197.    ((typep arg 'double-float)
  198.     (%random-double-float arg state))
  199.    ((integerp arg)
  200.     (%random-integer arg state))
  201.    (t
  202.     (error 'simple-type-error :expected-type '(real (0)) :datum arg
  203.        :format-control "Argument is not a positive real number: ~S"
  204.        :format-arguments (list arg)))))
  205.  
  206.  
  207.  
  208.  
  209.  
  210.  
  211.  
  212.  
  213.  
  214.  
  215.  
  216.